0. set up

# source the cleaning/summarizing script, which sources the set up
source(here::here("code", "01_cleaning-and-summarizing.R"))

1. Biomass

a. Guiding questions

What’s the percent community biomass that can be attributed to:
1. Laminariales, and
2. the 11 abundant species from Miller et al.?

Data question: what is the difference between the percent cover data and the percent cover that’s in the column in the biomass data?

Stray thoughts: introduced/native species trait column?

b. Test case: SCTW from 3 different sampling points

To make sure the code was working, I tested this out with “SCTW_2010-07-29”, “SCTW_2016-07-20”, and “SCTW_2020-07-29”. The issue comes from scenarios where species were counted at more than one transect - for example, at SCTW_2020-07-29, EC was at both transect 2 and 3, but I only want the dry/wet biomass for EC for the whole survey.

When I create the data frame, there are then “duplicate” observations for total species biomass for EC because it was present at both transects. To fix this, I filtered the data frame for unique observations of the total and percent biomass calculations and double checked the total dry percent to make sure everything added up to 1.

i. Calculations

This is just a test to make sure the calculations worked.

ii. Plots

These are tests of the plots.

c. Biomass data set

Since this worked, I felt comfortable doing this for the whole biomass data frame.

i. Calculations

Note: as of 2022-03-04, the code to create the prop_biomass data frame is in the 00_set-up.R script.

ii. Plots

Made some plot functions to keep things tidy.

Note: FIX THIS. YOU NEED TO REDO THESE PLOTS BECAUSE THE DATA FRAME THEY’RE BUILT ON DOESN’T TAKE THE AVERAGE AT THE SITE ANYMORE.

2. Tables of species contribution to biomass

This is a big table of percent species contribution to community biomass for each site across all sampling points. It’s useful for figuring out where to sample what. 05_species-selection-tables.R is a separate script with this code, essentially to play around with sites.

3. Urchin density at each site

Plots using stat_summary of mean urchin density at each site and a timeseries:

urchin_plot <- ggplot(urchins %>% filter(site %in% c("bull", "aque", "napl", "ivee", "mohk", "carp")), aes(x = reorder(site, density), y = density)) +
  stat_summary(aes(fill = site), 
               fun = "mean", geom = "col") +
  stat_summary(aes(group = site),
               fun.data = mean_se, geom = "errorbar", width = 0.5) +
  labs(x = "Site", y = "Mean urchin density (number/m2)",
       title = "Purple urchin density") +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 35)) +
  theme_bw() +
  theme(legend.position = "none")

urchin_plot

urchins_time_plot <- urchins %>% 
  filter(site %in% c("bull", "aque", "napl", "ivee", "mohk", "carp")) %>% 
  filter(year > 2015) %>% 
  ggplot(aes(x = year, y = density)) +
  # stat_summary(aes(col = site), 
  #              fun = "sum", geom = "point", size = 3) +
  # stat_summary(aes(col = site),
  #              fun = "sum", geom = "line") +
  stat_summary(aes(col = site),
               fun = "mean", geom = "point", size = 3) +
  stat_summary(aes(col = site),
               fun = "mean", geom = "line") +
  stat_summary(aes(col = site),
               fun.data = mean_se, geom = "errorbar", width = 0.5) +
  labs(x = "Year", y = "Sum urchin density (number/m2)",
       title = "Purple urchin density through time") +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(-0.1, 15)) +
  theme_bw()

urchins_time_plot

Plot with pairwise comparisons of means:

# analysis of variance
urchin_aov <- aov(na.omit(urchins$density) ~ urchins$site)

# Tukey HSD
urchin_tukey <- TukeyHSD(urchin_aov)

# compact letter display
urchin_cld <- multcompLetters4(urchin_aov, urchin_tukey)

# output of compact letter display in data frame
urchin_cld_df <- enframe(urchin_cld[[1]]$Letters)

# summary data frame for plotting
urchin_plot_summary <- urchins %>% 
  # filter(site %in% c("bull", "aque", "ivee", "mohk", "napl", "carp")) %>% 
  group_by(site) %>% 
  summarize(mean = mean(density, na.rm = TRUE),
            se = se(density)) %>% 
  left_join(., urchin_cld_df, by = c("site" = "name")) %>% 
  # mutate(value = fct_relevel(value, "a", "b", "bc", "c")) %>% 
  arrange(value)
  # mutate(site = fct_relevel(site, "napl", "bull", "carp", "aque", "ivee", "mohk"))

# new plot
urchin_plot_new <- ggplot(urchin_plot_summary, aes(x = site, y = mean)) +
  geom_col() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
  geom_text(aes(label = value, y = mean + se + 1)) + 
  labs(x = "Site", y = "Mean urchin density (number/m2)",
       title = "Purple urchin density") +
  scale_y_continuous(expand = c(0, 0)) +
  # coord_cartesian(ylim = c(0, 28)) +
  theme_bw()

urchin_plot_new

4. Sand cover

Plot using stat_summary:

sand_plot <- ggplot(substrate %>% filter(common_name == "Sand"), aes(x = site, y = percent_cover)) +
  stat_summary(aes(fill = site), 
               fun = "mean", geom = "col") +
  stat_summary(aes(group = site), 
               fun.data = mean_se, geom = "errorbar", width = 0.5) +
  # facet_wrap(~site) +
  labs(x = "Site", y = "Mean sand percent cover (%)",
       title = "Sand cover across sites") +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 90)) +
  theme_bw() +
  theme(legend.position = "none")

sand_plot

sand_time_plot <- ggplot(substrate %>% filter(common_name == "Sand"), aes(x = year, y = percent_cover)) +
  stat_summary(aes(col = site), 
               fun = "mean", geom = "point", size = 3) +
  stat_summary(aes(col = site),
               fun = "mean", geom = "line") +
  stat_summary(aes(col = site), 
               fun.data = mean_se, geom = "errorbar", width = 0.5) +
  # facet_wrap(~site) +
  labs(x = "Site", y = "Mean sand percent cover (%)",
       title = "Sand cover across sites through time") +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(-1, 105)) +
  theme_bw() 

sand_time_plot

In case saving plots is of interest:

Plot with pairwise comparisons of means:

# filter out sand
sand <- substrate %>% 
  filter(common_name == "Sand")
  # filter(site %in% c("bull", "aque", "ivee", "mohk", "napl", "carp"))

# analysis of variance
sand_aov <- aov(na.omit(sand$percent_cover) ~ sand$site)

# Tukey HSD
sand_tukey <- TukeyHSD(sand_aov)

# compact letter display
sand_cld <- multcompLetters4(sand_aov, sand_tukey)

# results into data frame
sand_cld_df <- enframe(sand_cld[[1]]$Letters)

# summary data frame
sand_plot_summary <- sand %>% 
  group_by(site) %>% 
  summarize(mean = mean(percent_cover, na.rm = TRUE),
            se = se(percent_cover)) %>% 
  left_join(., sand_cld_df, by = c("site" = "name"))
  # mutate(site = fct_relevel(site, "napl", "bull", "carp", "aque", "ivee", "mohk"))

# new plot
sand_plot_new <- ggplot(sand_plot_summary, aes(x = site, y = mean)) +
  geom_col() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
  geom_text(aes(label = value, y = mean + se + 1)) + 
  labs(x = "Site", y = "Mean sand percent cover (%)",
       title = "Sand cover across sites") +
  scale_y_continuous(expand = c(0, 0)) +
  # coord_cartesian(ylim = c(0, 24)) +
  theme_bw()

sand_plot_new

5. Kelp cover

Plot using stat_summary:

biomass %>% 
  filter(sp_code == "MAPY") %>% 
  filter(!(site %in% c("scdi", "sctw"))) %>% 
  # replace all -99999 values with 0
  mutate(wm_gm2 = replace(wm_gm2, wm_gm2 < 0, 0)) %>% 
  ggplot(aes(x = reorder(site, wm_gm2), y = wm_gm2)) +
  stat_summary(aes(fill = site), 
               fun = "mean", geom = "col") +
  stat_summary(aes(group = site),
               fun.data = mean_se, geom = "errorbar", width = 0.5) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 7000)) +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "Site", y = "Kelp biomass (wet g/m2)",
       title = "Kelp biomass across sites")

biomass %>% 
  filter(sp_code == "MAPY") %>% 
  filter(site == "carp") %>% 
  # filter(!(site %in% c("scdi", "sctw"))) %>% 
  # replace all -99999 values with 0
  mutate(wm_gm2 = replace(wm_gm2, wm_gm2 < 0, 0)) %>% 
  ggplot(aes(x = year, y = wm_gm2)) +
  stat_summary(aes(col = site), 
               fun = "mean", geom = "point", size = 3) +
  stat_summary(aes(col = site),
               fun = "mean", geom = "line") +
  stat_summary(aes(col = site), 
               fun.data = mean_se, geom = "errorbar", width = 0.5) +
  theme_bw() 

Plot with pairwise comparisons:

Note: does it make sense that mean kelp biomass is the same across all sites?

mapy <- biomass %>% 
  # filter(!(site %in% c("sctw", "scdi"))) %>% 
  filter(site %in% c("aque", "ivee", "mohk", "napl", "carp")) %>% 
  filter(sp_code == "MAPY") %>% 
  group_by(site, year) %>% 
  summarize(total = sum(dry_gm2)) 
## `summarise()` has grouped output by 'site'. You can override using the `.groups` argument.
mapy_aov <- aov(na.omit(mapy$total) ~ mapy$site)

mapy_tukey <- TukeyHSD(mapy_aov)

mapy_cld <- multcompLetters4(mapy_aov, mapy_tukey)

mapy_cld_df <- enframe(mapy_cld[[1]]$Letters)

mapy_plot_summary <- mapy %>% 
  group_by(site) %>% 
  summarize(mean = mean(total, na.rm = TRUE),
            se = se(total)) %>% 
  left_join(., mapy_cld_df, by = c("site" = "name")) 

mapy_plot_new <- ggplot(mapy_plot_summary, aes(x = site, y = mean)) +
  geom_col() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
  geom_text(aes(label = value, y = mean + se + 70)) + 
  labs(x = "Site", y = "Mean biomass (g dry/m2)",
       title = "MAPY biomass across sites") +
  # scale_y_continuous(expand = c(0, 0)) +
  # coord_cartesian(ylim = c(0, 1600)) +
  theme_bw()

mapy_plot_new

6. Understory biomass

macrobio <- biomass %>% 
  filter(site %in% c("aque", "ivee", "mohk", "napl", "carp")) %>% 
  filter(sp_code != "MAPY") %>% 
  group_by(site, year) %>% 
  summarize(total = sum(dry_gm2))
## `summarise()` has grouped output by 'site'. You can override using the `.groups` argument.
macrobio_aov <- aov(macrobio$total ~ macrobio$site)

macrobio_tukey <- TukeyHSD(macrobio_aov)

macrobio_cld <- multcompLetters4(macrobio_aov, macrobio_tukey)

macrobio_cld_df <- enframe(macrobio_cld[[1]]$Letters)

macrobio_plot_summary <- macrobio %>% 
  group_by(site) %>% 
  summarize(mean = mean(total, na.rm = TRUE),
            se = sd(total)/sqrt(length(total))) %>% 
  left_join(., macrobio_cld_df, by = c("site" = "name")) 
  # mutate(site = fct_relevel(site, "napl", "carp", "aque", "ivee", "mohk"))

macrobio_plot_new <- ggplot(macrobio_plot_summary, aes(x = site, y = mean)) +
  geom_col() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
  geom_text(aes(label = value, y = mean + se + 50)) + 
  labs(x = "Site", y = "Mean biomass (g dry/m2)",
       title = "Understory macroalgal biomass across sites") +
  scale_y_continuous(expand = c(0, 0)) +
  # coord_cartesian(ylim = c(0, 1600)) +
  theme_bw()

macrobio_plot_new

# all algae

all_macro <- biomass %>% 
  filter(site %in% c("aque", "ivee", "mohk", "napl", "carp")) %>% 
  group_by(site, year) %>% 
  summarize(total = sum(dry_gm2))
## `summarise()` has grouped output by 'site'. You can override using the `.groups` argument.
all_macro_aov <- aov(all_macro$total ~ all_macro$site)

all_macro_tukey <- TukeyHSD(all_macro_aov)

all_macro_cld <- multcompLetters4(all_macro_aov, all_macro_tukey)

all_macro_cld_df <- enframe(all_macro_cld[[1]]$Letters)

all_macro_plot_summary <- all_macro %>% 
  group_by(site) %>% 
  summarize(mean = mean(total, na.rm = TRUE),
            se = sd(total)/sqrt(length(total))) %>% 
  left_join(., all_macro_cld_df, by = c("site" = "name")) 
  # mutate(site = fct_relevel(site, "napl", "carp", "aque", "ivee", "mohk"))

all_macro_plot_new <- ggplot(all_macro_plot_summary, aes(x = site, y = mean)) +
  geom_col() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
  geom_text(aes(label = value, y = mean + se + 80)) + 
  labs(x = "Site", y = "Mean biomass (g dry/m2)",
       title = "Whole community biomass across sites") +
  scale_y_continuous(expand = c(0, 0)) +
  # coord_cartesian(ylim = c(0, 4500)) +
  theme_bw()

all_macro_plot_new

7. Community metrics

a. wrangling

Going to make a site by species matrix for each sampling date. There are a couple of anomalies: CARP_2000-09-08 and NAPL_2000-09-18 don’t have any observations, and AQUE_2019-07-23 is missing the sand cover data because sand cover was taken on 2019-07-24. Going to fix that later.

Note: as of 2022-03-04, the code to create community_matrix and community_metadata are in the 00_set-up.R script.

b. Species richness/diversity

i. Richness

sp_rich <- specnumber(community_matrix) %>% 
  enframe() %>% 
  left_join(., community_metadata, by = c("name" = "sample_ID")) %>% 
  rename(richness = value)

sp_rich_aov <- aov(sp_rich$richness ~ sp_rich$site)

sp_rich_tukey <- TukeyHSD(sp_rich_aov)

sp_rich_cld <- multcompLetters4(sp_rich_aov, sp_rich_tukey)

sp_rich_cld_df <- enframe(sp_rich_cld[[1]]$Letters)

sp_rich_plot_summary <- sp_rich %>% 
  group_by(site) %>% 
  summarize(mean = mean(richness, na.rm = TRUE),
            se = sd(richness)/sqrt(length(richness))) %>% 
  left_join(., sp_rich_cld_df, by = c("site" = "name")) %>% 
  mutate(site = fct_relevel(site, "napl", "carp", "aque", "ivee", "mohk"))

sp_rich_plot_new <- ggplot(sp_rich_plot_summary, aes(x = site, y = mean)) +
  geom_col() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
  geom_text(aes(label = value, y = mean + se + 1)) + 
  labs(x = "Site", y = "Mean species richness",
       title = "Understory macroalgal species richness across sites") +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 25)) +
  theme_bw()

sp_rich_plot_new

ii. Diversity

sp_div <- diversity(community_matrix, index = "shannon") %>% 
  enframe() %>% 
  left_join(., community_metadata, by = c("name" = "sample_ID")) %>% 
  rename(shandiv = value)

sp_div_aov <- aov(sp_div$shandiv ~ sp_div$site)

sp_div_tukey <- TukeyHSD(sp_div_aov)

sp_div_cld <- multcompLetters4(sp_div_aov, sp_div_tukey)

sp_div_cld_df <- enframe(sp_div_cld[[1]]$Letters)

sp_div_plot_summary <- sp_div %>% 
  group_by(site) %>% 
  summarize(mean = mean(shandiv, na.rm = TRUE),
            se = sd(shandiv)/sqrt(length(shandiv))) %>% 
  left_join(., sp_div_cld_df, by = c("site" = "name")) %>% 
  mutate(site = fct_relevel(site, "napl", "carp", "aque", "ivee", "mohk"))

sp_div_plot_new <- ggplot(sp_div_plot_summary, aes(x = site, y = mean)) +
  geom_col() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
  geom_text(aes(label = value, y = mean + se + 0.2)) + 
  labs(x = "Site", y = "Mean Shannon diversity",
       title = "Understory macroalgal Shannon diversity across sites") +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 2.75)) +
  theme_bw()

sp_div_plot_new

sp_div %>% 
  ggplot(aes(x = site, y = shandiv)) +
  geom_violin() +
  geom_point() +
  geom_text(data = sp_div_cld_df, aes(x = name, y = 3, label = value))

c. SIMPER analysis

Wanted to know which species significantly contributed to differences in community composition.

sim <- with(community_metadata, simper(community_matrix, site)) 
head(summary(sim))
## $carp_napl
##                    average           sd      ratio          ava          avb    cumsum
## PTCA          2.455694e-01 0.2556539579 0.96055407  8.694819937 8.264586e+01 0.2861201
## CC            8.088812e-02 0.1036097668 0.78069975  5.306398810 2.057173e+01 0.3803653
## DL            7.458567e-02 0.1265995199 0.58914656 10.447767857 1.363259e+01 0.4672672
## CYOS          4.770110e-02 0.0789746402 0.60400527  8.626758988 5.942733e+00 0.5228451
## CO            4.735428e-02 0.0928932663 0.50977087 14.833928571 6.656250e-01 0.5780190
## R             4.551648e-02 0.0492601324 0.92400233  2.605208333 1.165536e+01 0.6310516
## EC            4.407481e-02 0.0601931871 0.73222256 11.176785714 6.742411e+00 0.6824044
## POLA          3.749029e-02 0.0957137185 0.39169195  8.316964286 2.410714e-02 0.7260855
## GR            3.347978e-02 0.0866076882 0.38656827 10.279017857 8.537946e-01 0.7650937
## RAT           3.026258e-02 0.0548837476 0.55139412  1.575297619 4.980134e+00 0.8003536
## LAFA          2.548957e-02 0.0530884206 0.48013421  0.010628472 6.305970e+00 0.8300522
## LS            1.959313e-02 0.0410853960 0.47688803  2.070833333 1.560417e+00 0.8528807
## DP            1.532958e-02 0.0345931828 0.44313886  4.162500000 4.901786e-01 0.8707417
## TALE          1.341341e-02 0.0489963948 0.27376316  2.595535714 1.607143e-02 0.8863700
## FB            1.161466e-02 0.0245496910 0.47310826  1.654910714 4.738839e-01 0.8999026
## CP            8.888653e-03 0.0327503077 0.27140671  0.843750000 7.955357e-01 0.9102590
## BR            7.159473e-03 0.0127488668 0.56157720  0.295238095 1.559226e+00 0.9186007
## Nandersoniana 6.707379e-03 0.0227786354 0.29445920  0.894940476 1.414583e+00 0.9264157
## BO            6.499400e-03 0.0204183675 0.31831142  0.901785714 8.416667e-01 0.9339883
## SCCA          6.424783e-03 0.0191308503 0.33583365  0.552455357 6.026786e-01 0.9414741
## EGME          5.748187e-03 0.0398427455 0.14427185  0.027107639 1.855239e+00 0.9481714
## CAL           5.699457e-03 0.0128044920 0.44511390  0.120238095 1.127232e+00 0.9548120
## CRYP          4.006483e-03 0.0314669032 0.12732373  2.049702381 0.000000e+00 0.9594801
## BF            3.767757e-03 0.0207867151 0.18125794  0.635119048 2.886905e-01 0.9638700
## DU            3.638079e-03 0.0166719754 0.21821525  1.478571429 0.000000e+00 0.9681089
## BLD           3.435456e-03 0.0093860603 0.36601687  0.483333333 3.475397e-01 0.9721116
## AMZO          3.360020e-03 0.0144763468 0.23210417  1.112202381 0.000000e+00 0.9760265
## COF           3.099095e-03 0.0105681556 0.29324841  0.000000000 6.026786e-01 0.9796373
## ER            2.247823e-03 0.0056742276 0.39614611  0.387946429 1.448661e-01 0.9822563
## SAMU          1.902505e-03 0.0255191202 0.07455214  0.007921627 3.156746e-01 0.9844730
## CG            1.500490e-03 0.0048250891 0.31097671  0.029910714 3.290179e-01 0.9862213
## CF            1.421864e-03 0.0031479900 0.45167369  0.072916667 3.541667e-01 0.9878779
## UEC           1.301632e-03 0.0128091952 0.10161704  0.269419643 1.584821e-02 0.9893945
## GS            1.118705e-03 0.0055895448 0.20014235  0.100446429 1.674107e-01 0.9906979
## GYSP          8.693371e-04 0.0032248178 0.26957711  0.054464286 1.543155e-01 0.9917108
## UV            6.940577e-04 0.0030300964 0.22905465  0.144642857 1.607143e-02 0.9925195
## HAGL          6.369483e-04 0.0051652848 0.12331330  0.000000000 1.446429e-01 0.9932616
## PL            6.322902e-04 0.0037412891 0.16900330  0.138318452 6.287202e-02 0.9939983
## FR            5.986664e-04 0.0033116541 0.18077565  0.107142857 1.071429e-01 0.9946958
## PRSP          5.896792e-04 0.0029515596 0.19978563  0.088020833 8.802083e-02 0.9953829
## AU            5.667135e-04 0.0038155426 0.14852763  0.171875000 0.000000e+00 0.9960432
## BRA           5.341714e-04 0.0023441884 0.22787050  0.063541667 8.169643e-02 0.9966655
## STIN          5.236085e-04 0.0030633004 0.17092954  0.118005952 7.261905e-02 0.9972756
## CZ            4.396347e-04 0.0029117167 0.15098816  0.138318452 0.000000e+00 0.9977878
## UBB           3.847832e-04 0.0028032770 0.13726194  0.064583333 0.000000e+00 0.9982362
## LI            3.363351e-04 0.0028926203 0.11627351  0.110937500 1.584821e-02 0.9986280
## LX            2.976988e-04 0.0029800372 0.09989769  0.043750000 1.458333e-02 0.9989749
## BPSE          2.271476e-04 0.0021710678 0.10462482  0.000000000 4.821429e-02 0.9992396
## SAHO          1.373665e-04 0.0007833223 0.17536392  0.035117064 8.591268e-04 0.9993996
## FG            1.170322e-04 0.0010029952 0.11668269  0.011904762 1.190476e-02 0.9995360
## FTHR          1.082206e-04 0.0008847429 0.12231873  0.000000000 2.083333e-02 0.9996621
## SAFU          9.881098e-05 0.0008503173 0.11620483  0.031250000 5.208333e-03 0.9997772
## EA            9.535253e-05 0.0009729069 0.09800787  0.000000000 1.468304e-02 0.9998883
## IR            7.439040e-05 0.0010242725 0.07262754  0.000000000 1.257440e-02 0.9999750
## PHSE          2.149273e-05 0.0002862812 0.07507559  0.009077381 0.000000e+00 1.0000000
## FASP          0.000000e+00 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## NEO           0.000000e+00 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## SELO          0.000000e+00 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## 
## $carp_bull
##                    average           sd      ratio          ava          avb    cumsum
## PTCA          3.090448e-01 0.2289500398 1.34983527  8.694819937 135.07343809 0.3373353
## BF            9.197659e-02 0.0935287221 0.98340477  0.635119048  38.49206349 0.4377316
## CYOS          7.641745e-02 0.0858260143 0.89037628  8.626758988  24.61623939 0.5211444
## DL            5.730248e-02 0.0861525975 0.66512771 10.447767857  16.36428571 0.5836924
## CO            4.502185e-02 0.0635774250 0.70814207 14.833928571   8.28333333 0.6328356
## Nandersoniana 3.806181e-02 0.0390967152 0.97352958  0.894940476  15.31984127 0.6743817
## EC            3.673142e-02 0.0396696416 0.92593280 11.176785714   0.51071429 0.7144756
## EGME          3.428148e-02 0.0471466249 0.72712479  0.027107639  12.71089601 0.7518952
## CC            3.417775e-02 0.0385564501 0.88643409  5.306398810  10.56250000 0.7892017
## POLA          2.826186e-02 0.0632361257 0.44692592  8.316964286   2.25000000 0.8200507
## GS            2.662457e-02 0.0483545079 0.55061199  0.100446429   9.01785714 0.8491125
## GR            2.208290e-02 0.0626752735 0.35233834 10.279017857   0.26785714 0.8732169
## R             1.505678e-02 0.0173127522 0.86969305  2.605208333   4.35714286 0.8896520
## AU            9.361261e-03 0.0259640187 0.36054747  0.171875000   1.77777778 0.8998702
## DP            9.322010e-03 0.0237906499 0.39183502  4.162500000   0.00000000 0.9100456
## TALE          8.699694e-03 0.0303837829 0.28632689  2.595535714   0.19285714 0.9195416
## FB            6.980982e-03 0.0144868999 0.48188240  1.654910714   0.43869048 0.9271617
## BO            6.928598e-03 0.0125652895 0.55140776  0.901785714   2.24444444 0.9347245
## LS            6.914296e-03 0.0170950689 0.40446141  2.070833333   0.00000000 0.9422718
## CRYP          6.521753e-03 0.0269508422 0.24198697  2.049702381   1.38571429 0.9493905
## RAT           5.641507e-03 0.0075841288 0.74385696  1.575297619   1.10337302 0.9555485
## BR            4.059314e-03 0.0077867227 0.52131231  0.295238095   1.40238095 0.9599794
## CF            3.273523e-03 0.0053257990 0.61465387  0.072916667   1.44444444 0.9635526
## BLD           3.057056e-03 0.0073184387 0.41771977  0.483333333   0.77026455 0.9668895
## DU            2.733843e-03 0.0129779054 0.21065366  1.478571429   0.00000000 0.9698736
## UV            2.558582e-03 0.0068386079 0.37413782  0.144642857   0.57857143 0.9726664
## PHSE          2.372076e-03 0.0052753723 0.44965083  0.009077381   0.87142857 0.9752556
## AMZO          2.338950e-03 0.0104913894 0.22294000  1.112202381   0.00000000 0.9778087
## SCCA          2.247923e-03 0.0065271360 0.34439647  0.552455357   0.17857143 0.9802624
## PL            2.243585e-03 0.0047384041 0.47348967  0.138318452   0.63710317 0.9827113
## CP            2.092718e-03 0.0080118151 0.26120399  0.843750000   0.00000000 0.9849956
## FTHR          1.915062e-03 0.0035805577 0.53485030  0.000000000   0.75000000 0.9870860
## NEO           1.836746e-03 0.0064514187 0.28470421  0.000000000   0.62500000 0.9890909
## CAL           1.222971e-03 0.0049922565 0.24497362  0.120238095   0.24047619 0.9904258
## ER            1.160133e-03 0.0031800399 0.36481718  0.387946429   0.02619048 0.9916921
## BRA           1.114865e-03 0.0027286734 0.40857386  0.063541667   0.38730159 0.9929090
## STIN          9.685879e-04 0.0028830848 0.33595540  0.118005952   0.31468254 0.9939663
## UEC           8.021044e-04 0.0081596284 0.09830158  0.269419643   0.04226190 0.9948418
## LX            7.337854e-04 0.0028315173 0.25914919  0.043750000   0.23333333 0.9956428
## PRSP          6.526429e-04 0.0029732107 0.21950779  0.088020833   0.20119048 0.9963552
## LAFA          6.249676e-04 0.0021643085 0.28876088  0.010628472   0.16287632 0.9970373
## CZ            5.641034e-04 0.0026609633 0.21199219  0.138318452   0.06706349 0.9976531
## SAFU          5.601568e-04 0.0016520556 0.33906653  0.031250000   0.13888889 0.9982645
## GYSP          4.853155e-04 0.0013990578 0.34688738  0.054464286   0.14523810 0.9987943
## FR            3.065458e-04 0.0019088689 0.16059029  0.107142857   0.04761905 0.9991289
## LI            2.271339e-04 0.0021488350 0.10570097  0.110937500   0.00000000 0.9993768
## UBB           2.179603e-04 0.0015475531 0.14084192  0.064583333   0.00000000 0.9996147
## IR            1.151722e-04 0.0009342512 0.12327751  0.000000000   0.03353175 0.9997404
## CG            8.870658e-05 0.0004796712 0.18493204  0.029910714   0.01329365 0.9998373
## SAHO          8.669998e-05 0.0005131600 0.16895310  0.035117064   0.00000000 0.9999319
## FG            3.518608e-05 0.0004979940 0.07065563  0.011904762   0.00000000 0.9999703
## SAMU          2.721045e-05 0.0003949935 0.06888834  0.007921627   0.00000000 1.0000000
## BPSE          0.000000e+00 0.0000000000        NaN  0.000000000   0.00000000 1.0000000
## COF           0.000000e+00 0.0000000000        NaN  0.000000000   0.00000000 1.0000000
## EA            0.000000e+00 0.0000000000        NaN  0.000000000   0.00000000 1.0000000
## FASP          0.000000e+00 0.0000000000        NaN  0.000000000   0.00000000 1.0000000
## HAGL          0.000000e+00 0.0000000000        NaN  0.000000000   0.00000000 1.0000000
## SELO          0.000000e+00 0.0000000000        NaN  0.000000000   0.00000000 1.0000000
## 
## $carp_mohk
##                    average           sd      ratio          ava        avb    cumsum
## PTCA          1.629987e-01 0.2123082228 0.76774558  8.694819937 44.0066146 0.1798691
## CO            9.507108e-02 0.1112406399 0.85464345 14.833928571 13.4456250 0.2847801
## CYOS          9.466860e-02 0.1104015005 0.85749376  8.626758988 13.9713325 0.3892470
## EC            8.792094e-02 0.1173889353 0.74897125 11.176785714  0.4331250 0.4862678
## R             7.908317e-02 0.0763753306 1.03545443  2.605208333 17.1562500 0.5735361
## CC            6.670531e-02 0.0640267406 1.04183512  5.306398810 13.7840625 0.6471454
## GR            4.800132e-02 0.1003171248 0.47849582 10.279017857  2.3906250 0.7001149
## POLA          4.388456e-02 0.1093719281 0.40124151  8.316964286  0.1012500 0.7485415
## DL            4.053225e-02 0.1070211051 0.37873137 10.447767857  0.8156250 0.7932689
## RAT           2.072370e-02 0.0241639399 0.85762927  1.575297619  4.1456250 0.8161375
## DP            1.640438e-02 0.0368369271 0.44532423  4.162500000  0.4387500 0.8342397
## TALE          1.632516e-02 0.0564391536 0.28925237  2.595535714  0.2025000 0.8522546
## LS            1.526779e-02 0.0401103140 0.38064501  2.070833333  0.0000000 0.8691026
## FB            1.430644e-02 0.0337190026 0.42428425  1.654910714  0.3196875 0.8848897
## BLD           1.271000e-02 0.0277770116 0.45757261  0.483333333  2.2765000 0.8989152
## BO            1.139346e-02 0.0266472330 0.42756637  0.901785714  1.5150000 0.9114879
## EGME          1.078418e-02 0.0243393811 0.44307554  0.027107639  2.1278385 0.9233883
## BR            7.688069e-03 0.0178998054 0.42950575  0.295238095  1.3175000 0.9318721
## GS            5.697721e-03 0.0151997199 0.37485696  0.100446429  1.2656250 0.9381595
## CRYP          4.895620e-03 0.0334951620 0.14615903  2.049702381  0.1212500 0.9435618
## DU            4.324915e-03 0.0179240565 0.24129109  1.478571429  0.0675000 0.9483344
## SCCA          4.130474e-03 0.0133412827 0.30960093  0.552455357  0.1406250 0.9528924
## GYSP          3.949196e-03 0.0101914712 0.38750005  0.054464286  0.6100000 0.9572503
## AMZO          3.806647e-03 0.0161119745 0.23626201  1.112202381  0.0000000 0.9614509
## Nandersoniana 3.702192e-03 0.0153584594 0.24105232  0.894940476  0.4850000 0.9655363
## CP            3.662915e-03 0.0142207179 0.25757595  0.843750000  0.0000000 0.9695783
## AU            3.542542e-03 0.0120498855 0.29398967  0.171875000  0.5468750 0.9734875
## CZ            3.366317e-03 0.0101240971 0.33250538  0.138318452  0.4753125 0.9772023
## PRSP          2.803459e-03 0.0090408161 0.31008917  0.088020833  0.3696875 0.9802959
## ER            2.278642e-03 0.0069553852 0.32760833  0.387946429  0.0000000 0.9828104
## BF            1.917702e-03 0.0123580476 0.15517841  0.635119048  0.0000000 0.9849265
## BRA           1.724404e-03 0.0051040139 0.33785261  0.063541667  0.3050000 0.9868294
## LAFA          1.590329e-03 0.0053810939 0.29554016  0.010628472  0.2301229 0.9885844
## PL            1.539362e-03 0.0068096210 0.22605690  0.138318452  0.2112500 0.9902830
## CAL           1.452365e-03 0.0056392407 0.25754623  0.120238095  0.1262500 0.9918857
## UEC           1.311896e-03 0.0142724374 0.09191814  0.269419643  0.0000000 0.9933334
## FR            1.215001e-03 0.0047748041 0.25446099  0.107142857  0.3250000 0.9946742
## LI            1.040688e-03 0.0048743876 0.21350122  0.110937500  0.1331250 0.9958226
## UV            9.769952e-04 0.0037999680 0.25710618  0.144642857  0.0337500 0.9969007
## UBB           7.340540e-04 0.0035975763 0.20404126  0.064583333  0.0775000 0.9977107
## STIN          4.335038e-04 0.0024839272 0.17452356  0.118005952  0.0381250 0.9981891
## CF            4.089273e-04 0.0012163657 0.33618776  0.072916667  0.0656250 0.9986403
## PHSE          3.729136e-04 0.0024715066 0.15088513  0.009077381  0.0381250 0.9990518
## SAFU          2.418915e-04 0.0013145763 0.18400720  0.031250000  0.0437500 0.9993188
## CG            2.110499e-04 0.0009967199 0.21174442  0.029910714  0.0209375 0.9995517
## SAHO          1.532170e-04 0.0008932640 0.17152489  0.035117064  0.0000000 0.9997207
## LX            1.291713e-04 0.0012608312 0.10244928  0.043750000  0.0000000 0.9998633
## FG            6.468797e-05 0.0009041190 0.07154807  0.011904762  0.0000000 0.9999347
## SAMU          5.921697e-05 0.0008700064 0.06806499  0.007921627  0.0000000 1.0000000
## BPSE          0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## COF           0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## EA            0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## FASP          0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## FTHR          0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## HAGL          0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## IR            0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## NEO           0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## SELO          0.000000e+00 0.0000000000        NaN  0.000000000  0.0000000 1.0000000
## 
## $carp_ivee
##                    average           sd      ratio          ava          avb    cumsum
## DL            0.1312661303 0.1900446250 0.69071214 10.447767857 2.624500e+01 0.1470734
## PTCA          0.1197036199 0.1761482632 0.67956174  8.694819937 1.998016e+01 0.2811918
## CYOS          0.0865615525 0.1105684061 0.78287782  8.626758988 1.390637e+01 0.3781772
## CO            0.0820668440 0.1242198862 0.66065786 14.833928571 7.484583e+00 0.4701267
## EC            0.0752750791 0.0890583815 0.84523296 11.176785714 1.705000e+00 0.5544665
## POLA          0.0615202171 0.1206187433 0.51003862  8.316964286 4.230000e+00 0.6233950
## GR            0.0335516730 0.0962484522 0.34859442 10.279017857 3.125000e-02 0.6609870
## BF            0.0323666054 0.0709527129 0.45617150  0.635119048 5.496667e+00 0.6972513
## GS            0.0299702168 0.0615556399 0.48688011  0.100446429 6.625000e+00 0.7308305
## CC            0.0286887442 0.0535456148 0.53578139  5.306398810 3.966806e+00 0.7629740
## R             0.0237801211 0.0444181865 0.53536902  2.605208333 3.117778e+00 0.7896178
## Nandersoniana 0.0225936897 0.0484140311 0.46667648  0.894940476 4.149444e+00 0.8149322
## TALE          0.0180237829 0.0547605608 0.32913803  2.595535714 5.550000e-01 0.8351265
## FB            0.0177747700 0.0339181192 0.52404940  1.654910714 1.228333e+00 0.8550417
## DP            0.0158521388 0.0366420429 0.43262159  4.162500000 1.650000e-01 0.8728028
## LS            0.0154228949 0.0362885040 0.42500774  2.070833333 1.088889e-01 0.8900829
## RAT           0.0108499994 0.0170140440 0.63770844  1.575297619 7.351389e-01 0.9022395
## EGME          0.0099991100 0.0339286238 0.29471016  0.027107639 2.511308e+00 0.9134427
## BR            0.0076581706 0.0204210909 0.37501281  0.295238095 9.300000e-01 0.9220231
## BO            0.0076310186 0.0240526797 0.31726272  0.901785714 7.294444e-01 0.9305730
## SAMU          0.0061392105 0.0494033702 0.12426704  0.007921627 2.982759e+00 0.9374515
## AU            0.0058244785 0.0133882138 0.43504523  0.171875000 1.234722e+00 0.9439774
## CRYP          0.0055297062 0.0352598296 0.15682737  2.049702381 1.077778e-01 0.9501730
## PHSE          0.0048949395 0.0184521368 0.26527765  0.009077381 7.794444e-01 0.9556574
## SCCA          0.0048623369 0.0135937524 0.35768909  0.552455357 2.812500e-01 0.9611052
## BLD           0.0046276268 0.0126573174 0.36560881  0.483333333 3.888148e-01 0.9662901
## DU            0.0040411007 0.0180548791 0.22382319  1.478571429 0.000000e+00 0.9708179
## AMZO          0.0039123107 0.0158838362 0.24630767  1.112202381 2.805556e-02 0.9752013
## CP            0.0036636905 0.0137773343 0.26592158  0.843750000 0.000000e+00 0.9793062
## CF            0.0026185647 0.0057238008 0.45748704  0.072916667 4.277778e-01 0.9822401
## ER            0.0024485328 0.0063627871 0.38482079  0.387946429 5.958333e-02 0.9849835
## BRA           0.0018300085 0.0083293562 0.21970588  0.063541667 2.033333e-01 0.9870338
## UV            0.0015325063 0.0059160217 0.25904339  0.144642857 1.500000e-01 0.9887509
## FTHR          0.0015199490 0.0051650240 0.29427724  0.000000000 3.013889e-01 0.9904539
## UEC           0.0014030451 0.0140551331 0.09982439  0.269419643 2.958333e-02 0.9920259
## UBB           0.0013207516 0.0060098354 0.21976503  0.064583333 2.238889e-01 0.9935057
## GYSP          0.0010254566 0.0039650690 0.25862264  0.054464286 1.355556e-01 0.9946546
## CAL           0.0006641022 0.0041075407 0.16167879  0.120238095 0.000000e+00 0.9953987
## FR            0.0006278653 0.0033488718 0.18748563  0.107142857 6.666667e-02 0.9961022
## CZ            0.0004967955 0.0032102137 0.15475465  0.138318452 0.000000e+00 0.9966588
## SAFU          0.0004838790 0.0020137357 0.24028924  0.031250000 8.750000e-02 0.9972009
## LX            0.0004374759 0.0025369251 0.17244338  0.043750000 5.444444e-02 0.9976911
## LI            0.0004015853 0.0031774971 0.12638417  0.110937500 2.958333e-02 0.9981410
## PL            0.0003666459 0.0023285160 0.15745905  0.138318452 0.000000e+00 0.9985518
## PRSP          0.0003335074 0.0026474804 0.12597162  0.088020833 2.347222e-02 0.9989255
## LAFA          0.0003151463 0.0019141828 0.16463750  0.010628472 4.222917e-02 0.9992786
## STIN          0.0003001699 0.0023514742 0.12765179  0.118005952 0.000000e+00 0.9996149
## SAHO          0.0001533313 0.0008698014 0.17628314  0.035117064 2.222222e-04 0.9997867
## CG            0.0001254836 0.0007117883 0.17629349  0.029910714 9.305556e-03 0.9999273
## FG            0.0000648734 0.0008859051 0.07322839  0.011904762 0.000000e+00 1.0000000
## BPSE          0.0000000000 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## COF           0.0000000000 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## EA            0.0000000000 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## FASP          0.0000000000 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## HAGL          0.0000000000 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## IR            0.0000000000 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## NEO           0.0000000000 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## SELO          0.0000000000 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## 
## $carp_aque
##                    average           sd      ratio          ava          avb    cumsum
## PTCA          1.505860e-01 0.1974155792 0.76278691  8.694819937 3.010669e+01 0.1697838
## DL            1.487785e-01 0.2109229092 0.70536925 10.447767857 3.133125e+01 0.3375296
## EC            9.373259e-02 0.1318370641 0.71097302 11.176785714 2.414547e+00 0.4432119
## CYOS          7.683385e-02 0.1115759669 0.68862363  8.626758988 9.861216e+00 0.5298411
## CO            6.821917e-02 0.1113665372 0.61256436 14.833928571 3.213362e+00 0.6067573
## POLA          6.812189e-02 0.1443652138 0.47187189  8.316964286 4.922845e+00 0.6835638
## GR            3.566193e-02 0.1042796397 0.34198366 10.279017857 0.000000e+00 0.7237722
## R             3.239649e-02 0.0460650028 0.70327774  2.605208333 4.272629e+00 0.7602988
## CC            2.486131e-02 0.0525042771 0.47351023  5.306398810 2.185345e+00 0.7883296
## TALE          2.096720e-02 0.0646640638 0.32424806  2.595535714 7.913793e-01 0.8119699
## LS            1.855623e-02 0.0477778732 0.38838533  2.070833333 6.336207e-02 0.8328918
## DP            1.797162e-02 0.0411647354 0.43657796  4.162500000 1.512931e-01 0.8531545
## FB            1.778685e-02 0.0395719638 0.44948109  1.654910714 5.298491e-01 0.8732090
## RAT           1.310815e-02 0.0230978745 0.56750477  1.575297619 8.663793e-01 0.8879882
## AU            1.104456e-02 0.0285461994 0.38690112  0.171875000 1.689655e+00 0.9004408
## GS            9.445987e-03 0.0247688436 0.38136569  0.100446429 1.575970e+00 0.9110910
## BO            9.042056e-03 0.0267112618 0.33851100  0.901785714 7.400862e-01 0.9212858
## BR            6.343935e-03 0.0141916302 0.44701944  0.295238095 9.219828e-01 0.9284386
## Nandersoniana 5.999693e-03 0.0200053150 0.29990497  0.894940476 7.943966e-01 0.9352031
## CRYP          5.872691e-03 0.0362177102 0.16214971  2.049702381 1.254310e-01 0.9418245
## CP            5.463515e-03 0.0189751600 0.28792985  0.843750000 2.443966e-01 0.9479846
## SCCA          4.829829e-03 0.0153858522 0.31391362  0.552455357 1.454741e-01 0.9534301
## BF            4.710111e-03 0.0159120540 0.29600897  0.635119048 5.853448e-01 0.9587407
## DU            4.279185e-03 0.0189715645 0.22555786  1.478571429 0.000000e+00 0.9635654
## AMZO          4.217860e-03 0.0176383907 0.23912950  1.112202381 0.000000e+00 0.9683210
## BLD           3.594995e-03 0.0112036398 0.32087741  0.483333333 2.250000e-01 0.9723743
## ER            2.716476e-03 0.0082045793 0.33109263  0.387946429 2.489224e-02 0.9754371
## UV            2.512290e-03 0.0068189010 0.36843037  0.144642857 3.840517e-01 0.9782697
## EGME          1.845034e-03 0.0128303836 0.14380196  0.027107639 2.172195e-01 0.9803499
## BRA           1.681693e-03 0.0072122150 0.23317283  0.063541667 2.103448e-01 0.9822460
## UEC           1.657701e-03 0.0160432854 0.10332680  0.269419643 2.295259e-02 0.9841151
## SAMU          1.499332e-03 0.0128068709 0.11707249  0.007921627 1.950359e-01 0.9858055
## CF            1.467704e-03 0.0038851556 0.37777215  0.072916667 2.338362e-01 0.9874604
## GYSP          1.200532e-03 0.0049878660 0.24069055  0.054464286 1.709052e-01 0.9888139
## SELO          1.167799e-03 0.0090619374 0.12886856  0.000000000 1.396552e-01 0.9901306
## CAL           1.043564e-03 0.0054220778 0.19246564  0.120238095 4.353448e-02 0.9913072
## PHSE          9.191940e-04 0.0037901100 0.24252437  0.009077381 2.366379e-01 0.9923436
## UBB           8.769062e-04 0.0045678448 0.19197373  0.064583333 8.017241e-02 0.9933323
## STIN          7.196498e-04 0.0055606242 0.12941888  0.118005952 3.943966e-02 0.9941437
## LI            5.673895e-04 0.0040027860 0.14174865  0.110937500 4.590517e-02 0.9947834
## FASP          5.374982e-04 0.0033386566 0.16099236  0.000000000 1.051724e-01 0.9953894
## CZ            5.341173e-04 0.0035014174 0.15254316  0.138318452 0.000000e+00 0.9959917
## FR            5.130059e-04 0.0031828609 0.16117760  0.107142857 5.172414e-02 0.9965701
## LX            4.781302e-04 0.0026632471 0.17952904  0.043750000 6.336207e-02 0.9971092
## NEO           4.775126e-04 0.0031700373 0.15063312  0.000000000 1.454741e-01 0.9976475
## PRSP          4.326203e-04 0.0029773153 0.14530552  0.088020833 3.642241e-02 0.9981353
## PL            4.299539e-04 0.0024933571 0.17243975  0.138318452 1.821121e-02 0.9986201
## FTHR          3.774483e-04 0.0016700850 0.22600542  0.000000000 8.297414e-02 0.9990456
## LAFA          1.859816e-04 0.0015872062 0.11717544  0.010628472 1.303376e-02 0.9992553
## SAHO          1.732297e-04 0.0010044809 0.17245696  0.035117064 8.620690e-05 0.9994507
## IR            1.413228e-04 0.0016593587 0.08516712  0.000000000 1.821121e-02 0.9996100
## CG            1.088585e-04 0.0007335932 0.14839084  0.029910714 0.000000e+00 0.9997327
## FG            1.084103e-04 0.0010890267 0.09954784  0.011904762 8.620690e-03 0.9998550
## SAFU          9.846824e-05 0.0009588661 0.10269237  0.031250000 0.000000e+00 0.9999660
## EA            3.016904e-05 0.0003365065 0.08965367  0.000000000 7.212644e-03 1.0000000
## BPSE          0.000000e+00 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## COF           0.000000e+00 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## HAGL          0.000000e+00 0.0000000000        NaN  0.000000000 0.000000e+00 1.0000000
## 
## $napl_bull
##                    average           sd      ratio          ava          avb    cumsum
## PTCA          2.776579e-01 2.080389e-01 1.33464424 8.264586e+01 135.07343809 0.3617611
## BF            8.074493e-02 8.432750e-02 0.95751603 2.886905e-01  38.49206349 0.4669639
## CYOS          5.971607e-02 7.076090e-02 0.84391334 5.942733e+00  24.61623939 0.5447680
## DL            5.356836e-02 7.445622e-02 0.71946116 1.363259e+01  16.36428571 0.6145624
## CC            4.729296e-02 5.755267e-02 0.82173349 2.057173e+01  10.56250000 0.6761805
## Nandersoniana 3.318287e-02 3.419358e-02 0.97044149 1.414583e+00  15.31984127 0.7194145
## EGME          3.210429e-02 4.576318e-02 0.70153107 1.855239e+00  12.71089601 0.7612433
## R             2.422633e-02 2.789667e-02 0.86843085 1.165536e+01   4.35714286 0.7928078
## GS            2.277136e-02 4.222268e-02 0.53931592 1.674107e-01   9.01785714 0.8224767
## CO            2.027309e-02 2.699765e-02 0.75092060 6.656250e-01   8.28333333 0.8488905
## EC            1.705371e-02 1.890369e-02 0.90213624 6.742411e+00   0.51071429 0.8711098
## LAFA          1.436343e-02 2.842505e-02 0.50530894 6.305970e+00   0.16287632 0.8898240
## RAT           1.319429e-02 2.137444e-02 0.61729274 4.980134e+00   1.10337302 0.9070149
## AU            7.111229e-03 1.983495e-02 0.35852019 0.000000e+00   1.77777778 0.9162801
## POLA          6.086456e-03 2.467037e-02 0.24671121 2.410714e-02   2.25000000 0.9242102
## BR            5.410701e-03 7.936182e-03 0.68177624 1.559226e+00   1.40238095 0.9312598
## BO            5.390379e-03 7.358635e-03 0.73252434 8.416667e-01   2.24444444 0.9382829
## LS            4.775916e-03 1.503448e-02 0.31766418 1.560417e+00   0.00000000 0.9445055
## CAL           3.269553e-03 7.043646e-03 0.46418477 1.127232e+00   0.24047619 0.9487654
## CF            3.112419e-03 4.764849e-03 0.65320408 3.541667e-01   1.44444444 0.9528205
## CRYP          3.107866e-03 8.153521e-03 0.38116855 0.000000e+00   1.38571429 0.9568698
## FB            3.038833e-03 9.027028e-03 0.33663711 4.738839e-01   0.43869048 0.9608291
## CP            2.635013e-03 1.416004e-02 0.18608804 7.955357e-01   0.00000000 0.9642623
## GR            2.506314e-03 6.053939e-03 0.41399716 8.537946e-01   0.26785714 0.9675277
## BLD           2.356245e-03 5.192892e-03 0.45374423 3.475397e-01   0.77026455 0.9705977
## SCCA          2.227789e-03 7.582121e-03 0.29382132 6.026786e-01   0.17857143 0.9735003
## PHSE          2.040969e-03 4.585201e-03 0.44512093 0.000000e+00   0.87142857 0.9761595
## UV            1.867275e-03 5.411747e-03 0.34504100 1.607143e-02   0.57857143 0.9785924
## PL            1.853173e-03 4.108275e-03 0.45108307 6.287202e-02   0.63710317 0.9810069
## FTHR          1.686795e-03 3.162811e-03 0.53332164 2.083333e-02   0.75000000 0.9832046
## NEO           1.565331e-03 5.568965e-03 0.28108120 0.000000e+00   0.62500000 0.9852441
## COF           1.533881e-03 5.049379e-03 0.30377609 6.026786e-01   0.00000000 0.9872426
## DP            1.393048e-03 5.969179e-03 0.23337351 4.901786e-01   0.00000000 0.9890576
## BRA           9.744933e-04 2.345415e-03 0.41548855 8.169643e-02   0.38730159 0.9903272
## SAMU          9.207575e-04 1.299543e-02 0.07085243 3.156746e-01   0.00000000 0.9915269
## TALE          8.511697e-04 4.662908e-03 0.18254054 1.607143e-02   0.19285714 0.9926359
## STIN          8.238719e-04 2.514940e-03 0.32759101 7.261905e-02   0.31468254 0.9937093
## CG            8.067520e-04 2.563852e-03 0.31466398 3.290179e-01   0.01329365 0.9947604
## GYSP          6.895111e-04 1.830906e-03 0.37659561 1.543155e-01   0.14523810 0.9956588
## LX            6.181498e-04 2.548389e-03 0.24256496 1.458333e-02   0.23333333 0.9964642
## PRSP          6.033651e-04 2.373910e-03 0.25416514 8.802083e-02   0.20119048 0.9972503
## SAFU          4.267935e-04 1.321708e-03 0.32291070 5.208333e-03   0.13888889 0.9978064
## ER            3.595011e-04 9.671849e-04 0.37169846 1.448661e-01   0.02619048 0.9982748
## HAGL          3.488862e-04 2.849984e-03 0.12241691 1.446429e-01   0.00000000 0.9987293
## FR            2.792231e-04 1.556732e-03 0.17936493 1.071429e-01   0.04761905 0.9990931
## CZ            2.122696e-04 1.437574e-03 0.14765819 0.000000e+00   0.06706349 0.9993697
## IR            1.332323e-04 9.488075e-04 0.14042080 1.257440e-02   0.03353175 0.9995433
## UEC           1.267877e-04 1.072501e-03 0.11821686 1.584821e-02   0.04226190 0.9997085
## BPSE          1.245917e-04 1.214452e-03 0.10259088 4.821429e-02   0.00000000 0.9998708
## EA            4.426948e-05 4.468248e-04 0.09907569 1.468304e-02   0.00000000 0.9999285
## FG            3.185882e-05 3.203910e-04 0.09943731 1.190476e-02   0.00000000 0.9999700
## LI            2.051270e-05 2.690408e-04 0.07624384 1.584821e-02   0.00000000 0.9999967
## SAHO          2.515191e-06 3.552632e-05 0.07079796 8.591268e-04   0.00000000 1.0000000
## AMZO          0.000000e+00 0.000000e+00        NaN 0.000000e+00   0.00000000 1.0000000
## DU            0.000000e+00 0.000000e+00        NaN 0.000000e+00   0.00000000 1.0000000
## FASP          0.000000e+00 0.000000e+00        NaN 0.000000e+00   0.00000000 1.0000000
## SELO          0.000000e+00 0.000000e+00        NaN 0.000000e+00   0.00000000 1.0000000
## UBB           0.000000e+00 0.000000e+00        NaN 0.000000e+00   0.00000000 1.0000000

d. NMDS

i. Jaccard (presence/absence)

community_mds_jacc <- metaMDS(community_presabs, distance = "jaccard")
## Run 0 stress 0.2324413 
## Run 1 stress 0.2333624 
## Run 2 stress 0.2500419 
## Run 3 stress 0.2528057 
## Run 4 stress 0.2347918 
## Run 5 stress 0.2339213 
## Run 6 stress 0.2327236 
## ... Procrustes: rmse 0.006726982  max resid 0.09394191 
## Run 7 stress 0.2350828 
## Run 8 stress 0.2372256 
## Run 9 stress 0.2339764 
## Run 10 stress 0.2334229 
## Run 11 stress 0.2368917 
## Run 12 stress 0.2389732 
## Run 13 stress 0.2364392 
## Run 14 stress 0.2327612 
## ... Procrustes: rmse 0.006176838  max resid 0.07131762 
## Run 15 stress 0.2418717 
## Run 16 stress 0.240909 
## Run 17 stress 0.2332244 
## Run 18 stress 0.2495368 
## Run 19 stress 0.2332678 
## Run 20 stress 0.2381245 
## *** No convergence -- monoMDS stopping criteria:
##      4: no. of iterations >= maxit
##     13: stress ratio > sratmax
##      3: scale factor of the gradient < sfgrmin
community_dist_jacc <- vegdist(community_presabs, method = "jaccard")

stressplot(community_mds_jacc)

community_mds_sites_jacc <- scores(community_mds_jacc, display = "sites") %>% 
  as.data.frame() %>% 
  rownames_to_column("sample_ID") %>% 
  left_join(., community_metadata, by = "sample_ID") 

community_mds_spp_jacc <- scores(community_mds_jacc, display = "species") %>% 
  as.data.frame() %>% 
  rownames_to_column("sp_code") %>% 
  left_join(., coarse_traits, by = "sp_code") %>% 
  filter(sp_code %in% algae_proposal)

perm <- adonis2(community_dist_jacc ~ site, data = community_metadata)

perm
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = community_dist_jacc ~ site, data = community_metadata)
##           Df SumOfSqs      R2      F Pr(>F)    
## site       5   29.442 0.15199 22.906  0.001 ***
## Residual 639  164.265 0.84801                  
## Total    644  193.707 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
disper <- betadisper(community_dist_jacc, community_metadata$site)

anova(disper)
## Analysis of Variance Table
## 
## Response: Distances
##            Df Sum Sq  Mean Sq F value    Pr(>F)    
## Groups      5 1.2762 0.255236   32.14 < 2.2e-16 ***
## Residuals 639 5.0745 0.007941                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(disper)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = distances ~ group, data = df)
## 
## $group
##                   diff          lwr          upr     p adj
## bull-aque -0.106271846 -0.146136993 -0.066406699 0.0000000
## carp-aque  0.015459985 -0.015289770  0.046209740 0.7042879
## ivee-aque  0.012470074 -0.023310680  0.048250828 0.9191448
## mohk-aque -0.053008215 -0.099713866 -0.006302564 0.0156226
## napl-aque -0.071037639 -0.101787394 -0.040287884 0.0000000
## carp-bull  0.121731831  0.084100706  0.159362956 0.0000000
## ivee-bull  0.118741920  0.076899161  0.160584679 0.0000000
## mohk-bull  0.053263631  0.001766330  0.104760932 0.0378176
## napl-bull  0.035234207 -0.002396918  0.072865332 0.0814327
## ivee-carp -0.002989911 -0.036263533  0.030283711 0.9998482
## mohk-carp -0.068468200 -0.113282143 -0.023654257 0.0002122
## napl-carp -0.086497624 -0.114290051 -0.058705196 0.0000000
## mohk-ivee -0.065478289 -0.113882879 -0.017073698 0.0016896
## napl-ivee -0.083507713 -0.116781335 -0.050234091 0.0000000
## napl-mohk -0.018029424 -0.062843367  0.026784519 0.8601005
fit_jacc <- community_metadata %>% 
  select(sample_ID, site, urchin_density, mean_sand, mean_daily_umol, kelp_dry) %>% 
  mutate(across(.cols = 3:4, ~replace(., is.na(.), 0))) %>% 
  envfit(community_mds_jacc, ., perm = 999, na.rm = TRUE)

fit_vect_jacc <- scores(fit, display = "vectors") %>% 
  as.data.frame()

nmds_plot_jacc <- ggplot() +
  # x and y axes
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  # points for sites
  geom_point(data = community_mds_sites_jacc, aes(x = NMDS1, y = NMDS2, color = site)) +
  stat_ellipse(data = community_mds_sites_jacc, aes(x = NMDS1, y = NMDS2, color = site), level = 0.95, size = 1) +
  # # arrows for species
  # geom_segment(data = community_mds_spp_jacc, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.2, "cm")), color = "blue") +
  # geom_text(data = community_mds_spp_jacc, aes(x = NMDS1, y = NMDS2, label = scientific_name), color = "blue") +
  # # arrows for environmental characteristics (vector fitting)
  geom_segment(data = fit_vect_jacc, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.2, "cm"))) +
  geom_text(data = fit_vect_jacc, aes(x = NMDS1, y = NMDS2, label = rownames(fit_vect_jacc))) +
  # theme
  theme_bw() +
  guides(colour = guide_legend(nrow = 5)) +
  theme(legend.position = "bottom", legend.box = "vertical", legend.margin = margin()) 

nmds_plot_jacc

ii. Bray-Curtis

community_mds_bray <- metaMDS(community_matrix, distance = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2397652 
## Run 1 stress 0.2460827 
## Run 2 stress 0.2590958 
## Run 3 stress 0.243352 
## Run 4 stress 0.2495171 
## Run 5 stress 0.2458121 
## Run 6 stress 0.2492311 
## Run 7 stress 0.2471866 
## Run 8 stress 0.2478991 
## Run 9 stress 0.24517 
## Run 10 stress 0.2592116 
## Run 11 stress 0.2403553 
## Run 12 stress 0.2508534 
## Run 13 stress 0.2502277 
## Run 14 stress 0.2489769 
## Run 15 stress 0.249233 
## Run 16 stress 0.2485628 
## Run 17 stress 0.2489904 
## Run 18 stress 0.2502947 
## Run 19 stress 0.2506252 
## Run 20 stress 0.247575 
## *** No convergence -- monoMDS stopping criteria:
##      4: no. of iterations >= maxit
##     16: stress ratio > sratmax
community_dist_bray <- vegdist(community_matrix, method = "bray")

stressplot(community_mds_bray)

community_mds_sites_bray <- scores(community_mds_bray, display = "sites") %>% 
  as.data.frame() %>% 
  rownames_to_column("sample_ID") %>% 
  left_join(., community_metadata, by = "sample_ID") 

community_mds_spp_bray <- scores(community_mds_bray, display = "species") %>% 
  as.data.frame() %>% 
  rownames_to_column("sp_code") %>% 
  left_join(., coarse_traits, by = "sp_code") %>% 
  filter(sp_code %in% algae_proposal)

perm <- adonis2(community_dist_bray ~ site, data = community_metadata)

perm
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = community_dist_bray ~ site, data = community_metadata)
##           Df SumOfSqs      R2      F Pr(>F)    
## site       5   36.368 0.15729 23.853  0.001 ***
## Residual 639  194.856 0.84271                  
## Total    644  231.225 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
disper <- betadisper(community_dist_bray, community_metadata$site)

anova(disper)
## Analysis of Variance Table
## 
## Response: Distances
##            Df Sum Sq  Mean Sq F value    Pr(>F)    
## Groups      5 1.4579 0.291571  25.107 < 2.2e-16 ***
## Residuals 639 7.4208 0.011613                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(disper)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = distances ~ group, data = df)
## 
## $group
##                   diff         lwr          upr     p adj
## bull-aque -0.166971431 -0.21517981 -0.118763050 0.0000000
## carp-aque -0.028129109 -0.06531437  0.009056153 0.2570883
## ivee-aque -0.007283007 -0.05055219  0.035986173 0.9968071
## mohk-aque -0.077425080 -0.13390559 -0.020944570 0.0013792
## napl-aque -0.065291455 -0.10247672 -0.028106194 0.0000100
## carp-bull  0.138842322  0.09333551  0.184349131 0.0000000
## ivee-bull  0.159688425  0.10908854  0.210288305 0.0000000
## mohk-bull  0.089546351  0.02727136  0.151821338 0.0006349
## napl-bull  0.101679976  0.05617317  0.147186784 0.0000000
## ivee-carp  0.020846102 -0.01939124  0.061083442 0.6766880
## mohk-carp -0.049295972 -0.10348886  0.004896921 0.0986720
## napl-carp -0.037162347 -0.07077135 -0.003553341 0.0203913
## mohk-ivee -0.070142074 -0.12867709 -0.011607059 0.0085374
## napl-ivee -0.058008449 -0.09824579 -0.017771109 0.0006074
## napl-mohk  0.012133625 -0.04205927  0.066326518 0.9879368
fit <- community_metadata %>% 
  select(sample_ID, site, urchin_density, mean_sand, mean_daily_umol, kelp_dry) %>% 
  mutate(across(.cols = 3:4, ~replace(., is.na(.), 0))) %>% 
  envfit(community_mds_bray, ., perm = 999, na.rm = TRUE)

fit_vect_bray <- scores(fit, display = "vectors") %>% 
  as.data.frame()

nmds_plot_bray <- ggplot() +
  geom_point(data = community_mds_sites_bray, aes(x = NMDS1, y = NMDS2, color = site, shape = site)) +
  # scale_color_manual(values = pal) +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  stat_ellipse(data = community_mds_sites_bray, aes(x = NMDS1, y = NMDS2, color = site), level = 0.95, size = 1) +
  geom_segment(data = community_mds_spp_bray, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.2, "cm"))) +
  geom_text(data = community_mds_spp_bray, aes(x = NMDS1, y = NMDS2, label = scientific_name)) +
  # geom_segment(data = fit_vect_bray, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.2, "cm"))) +
  # geom_text(data = fit_vect_bray, aes(x = NMDS1, y = NMDS2, label = rownames(fit_vect_bray))) +
  theme_bw() +
  guides(colour = guide_legend(nrow = 5)) +
  theme(legend.position = "bottom", legend.box = "vertical", legend.margin = margin()) 

nmds_plot_bray

8. Trait calculations

a. proportion of traits

“Mean proportion of biomass” for each species per trait across sampling dates: matrix where row is sampling date, and columns are traits spread across, and cells are total biomass that can be “assigned” to that trait. The prop_biomass data frame is already set up.

i. growth form

Going to test this out with a couple of sites just to make sure the math is right.

Now that it seems ok, going to do the rest of the data set.

prop_traits <- prop_biomass %>% 
  filter(site %in% c("aque", "napl", "ivee", "mohk", "carp")) %>%
  select(sample_ID, site, year, sp_code, growth_form,
         dry_gm2, wm_gm2, total_dry, total_wet, percent_sp_dry, percent_sp_wet) %>% 
  # create a new column where total biomass for the whole survey is calculated
  # group_by(site, year) %>% 
  # mutate(total_dry = sum(dry_gm2),
  #        total_wet = sum(wm_gm2)) %>% 
  # ungroup() %>% 
  # calculate biomass within each growth form
  group_by(sample_ID, growth_form) %>% 
  mutate(total_gf_dry = sum(dry_gm2),
         total_gf_wet = sum(wm_gm2)) %>% 
  ungroup() %>% 
  select(sample_ID, growth_form, total_gf_dry, total_dry, total_gf_wet, total_wet) %>% 
  unique() %>% 
  group_by(sample_ID, growth_form) %>% 
  summarize(percent_gf_dry = total_gf_dry/total_dry,
         percent_gf_wet = total_gf_wet/total_wet) %>% 
  ungroup() %>% 
  # double check: should sum to 1
  # group_by(sample_ID) %>%
  # mutate(test = sum(unique(percent_gf_dry)))
  # make the data frame wider
  select(sample_ID, growth_form, percent_gf_dry) %>% 
  unique() %>% 
  pivot_wider(names_from = "growth_form", values_from = "percent_gf_dry")
## `summarise()` has grouped output by 'sample_ID'. You can override using the `.groups` argument.

GLM???

9. Community-weighted means

These are calculated using the functcomp function in the FD package. The function requires a matrix with functional traits and a matrix of species abundances (or presence-absence). In the trait matrix, rows are species and columns are traits. In the species matrix, rows are samples and columns are species.

The first step is to get the coarse_traits data frame into a functional form (lol).

trait_matrix <- coarse_traits %>% 
  select(sp_code, growth_morph, growth_form, pigment_type, 
         height_class, longevity, posture, branching_yn) %>% 
  filter(sp_code != "MAPY") %>% 
  column_to_rownames("sp_code")

# making sure the rownames in the trait matrix == columns in the community matrix
rownames(trait_matrix) == colnames(community_matrix)
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [22] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Then the community matrix has to be a matrix, not a data frame (why?).

comm_matrix <- as.matrix(community_matrix)

Then, you’re ready to calculate community-weighted means.

CWM_cat <- functcomp(trait_matrix, comm_matrix, 
                  CWM.type = c("all")) %>% 
  rownames_to_column("sample_ID") %>% 
  left_join(., community_metadata, by = "sample_ID") %>% 
  # to get the glmm to work!!!!!!
  mutate(across(.cols = everything(), .fns = ~ replace(., . == 1, 0.999)))

ggplot(CWM_cat, aes(x = year, y = mean_daily_umol)) +
  geom_point(aes(col = site)) +
  geom_smooth(method = "lm", aes(col = site))
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 262 rows containing non-finite values (stat_smooth).
## Warning: Removed 262 rows containing missing values (geom_point).

For categorical traits, the output is the exact same as the prop_traits data frame up top - each cell is the proportion of the community biomass that is attributed to a specific category (for example, within growth_morph, at CARP in 2000, there is 92% solitary and 8 % aggregated).

Now to plot?

a. GLMs

Kill me!

lm_plot <- ggplot(CWM_cat, 
                  aes(x = mean_sand, y = growth_form_leathery_macrophyte)) +
  # scale_x_continuous(expand = c(0, 0)) +
  # scale_y_continuous(expand = c(0, 0)) +
  # coord_cartesian(ylim = c(-0.001, 1)) +
  geom_point(aes(col = site)) +
  geom_smooth(method = "lm", se = FALSE, aes(col = site)) +
  labs(x = "Mean sand cover", y = "Proportion biomass (dry g/m2)",
       title = "Proportion leathery macrophytes")
# missing values for mean sand

lm_plot
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing missing values (geom_point).

ggplot(CWM_cat, aes(x = growth_form_leathery_macrophyte)) +
  geom_histogram(binwidth = 0.1)

model_lm_01a <- glmmTMB(growth_form_leathery_macrophyte ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year), 
              data = CWM_cat, ziformula = ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry, family = beta_family(link = "logit"))
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN function evaluation
model_lm_01b <- glmmTMB(growth_form_leathery_macrophyte ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|year), 
              data = CWM_cat, ziformula = ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry, family = beta_family(link = "logit"))
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN function evaluation
model_lm_01c <- glmmTMB(growth_form_leathery_macrophyte ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + site + (1|year), 
              data = CWM_cat, ziformula = ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry, family = beta_family(link = "logit"))
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN function evaluation
# plot residuals
model_lm_01c_resid <- simulateResiduals(model_lm_01c)
plot(model_lm_01c_resid)

model_lm_nozi_resid <- simulateResiduals(model_lm_nozi)
plot(model_lm_nozi_resid)

# look at the model
model_lm_01c
## Formula:          growth_form_leathery_macrophyte ~ urchin_density + mean_sand +  
##     mean_daily_umol + kelp_dry + site + (1 | year)
## Zero inflation:                                   ~urchin_density + mean_sand + mean_daily_umol + kelp_dry
## Data: CWM_cat
##       AIC       BIC    logLik  df.resid 
##  50.04870 113.09144  -9.02435       364 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups Name        Std.Dev.
##  year   (Intercept) 0.2156  
## 
## Number of obs: 380 / Conditional model: year, 14
## 
## Dispersion parameter for beta family (): 3.17 
## 
## Fixed Effects:
## 
## Conditional model:
##     (Intercept)   urchin_density        mean_sand  mean_daily_umol         kelp_dry         sitecarp  
##      -1.096e+00       -9.588e-04        1.950e-02       -2.229e-06        5.067e-04       -4.798e-01  
##        siteivee         sitemohk         sitenapl  
##      -1.938e-01        1.267e+00        1.425e+00  
## 
## Zero-inflation model:
##     (Intercept)   urchin_density        mean_sand  mean_daily_umol         kelp_dry  
##      -2.309e+00        9.233e-02       -1.312e-01       -4.245e-05       -1.509e-03
summary(model_lm_01c)
##  Family: beta  ( logit )
## Formula:          growth_form_leathery_macrophyte ~ urchin_density + mean_sand +  
##     mean_daily_umol + kelp_dry + site + (1 | year)
## Zero inflation:                                   ~urchin_density + mean_sand + mean_daily_umol + kelp_dry
## Data: CWM_cat
## 
##      AIC      BIC   logLik deviance df.resid 
##     50.0    113.1     -9.0     18.0      364 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  year   (Intercept) 0.04646  0.2156  
## Number of obs: 380, groups:  year, 14
## 
## Dispersion parameter for beta family (): 3.17 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.096e+00  2.119e-01  -5.175 2.28e-07 ***
## urchin_density  -9.588e-04  7.810e-03  -0.123 0.902296    
## mean_sand        1.950e-02  4.554e-03   4.282 1.85e-05 ***
## mean_daily_umol -2.229e-06  1.298e-05  -0.172 0.863629    
## kelp_dry         5.067e-04  1.506e-04   3.365 0.000766 ***
## sitecarp        -4.798e-01  2.173e-01  -2.208 0.027230 *  
## siteivee        -1.938e-01  2.142e-01  -0.905 0.365540    
## sitemohk         1.267e+00  2.405e-01   5.270 1.37e-07 ***
## sitenapl         1.425e+00  2.143e-01   6.650 2.92e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.309e+00  3.801e-01  -6.073 1.25e-09 ***
## urchin_density   9.233e-02  1.577e-02   5.855 4.76e-09 ***
## mean_sand       -1.312e-01  9.638e-02  -1.361    0.174    
## mean_daily_umol -4.245e-05  3.670e-05  -1.157    0.247    
## kelp_dry        -1.509e-03  7.035e-04  -2.145    0.032 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# significant effects of sand and kelp biomass, not of urchin density and year
# why is it an increasing slope? this doesn't make sense

AIC(model_lm_01a, model_lm_01b, model_lm_01c)
##              df       AIC
## model_lm_01a 13  66.42352
## model_lm_01b 12 188.27592
## model_lm_01c 16  50.04870
# plot predicted values on top of predictor
temp <- ggpredict(model_lm, terms = "urchin_density")

# predict(model_lm)
plot(temp, add.data = TRUE)

b. Checking correlations

Checking to see if urchins, sand, kelp, and irradiance are at all correlated.

ggplot(CWM_cat, aes(x = urchin_density, y = mean_sand)) +
  geom_point(aes(col = site)) +
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing missing values (geom_point).

# create correlation matrix
cor_mat <- cor(CWM_cat[, c("urchin_density", "mean_sand", "mean_daily_umol", "kelp_dry")], use = "na.or.complete")
# plot correlation matrix
corrplot(cor_mat, type = "lower", addCoef.col = "black")

# 1) sand and urchins are slightly negatively correlated, 
# 2) irradiance is moderately positively correlated with urchin density, 
# 3) urchin density is moderately negatively correlated with kelp biomass, 
# 4) sand cover is slightly negatively correlated with kelp biomass, and 
# 5) average daily irradiance is slightly negatively correlated with kelp

# check for collinearity in the model
check_collinearity(model_lm_01c)
## # Check for Multicollinearity
## 
## * conditional component:
## 
## Low Correlation
## 
##             Term  VIF Increased SE Tolerance
##   urchin_density 1.34         1.16      0.74
##        mean_sand 1.53         1.24      0.65
##  mean_daily_umol 1.06         1.03      0.95
##         kelp_dry 1.62         1.27      0.62
##             site 2.27         1.51      0.44
## 
## * zero inflated component:
## 
## Low Correlation
## 
##             Term  VIF Increased SE Tolerance
##   urchin_density 1.16         1.08      0.86
##        mean_sand 1.02         1.01      0.98
##  mean_daily_umol 1.03         1.01      0.97
##         kelp_dry 1.12         1.06      0.89
# plot VIF for each variable (if >5, inflates)
plot(check_collinearity(model_lm_01c))

10. CWM for categorical and continuous traits

a. fixed CWM using FD::dbFD()

Use community_matrix and tbspp_matrix from set up script.

cwm_commat <- community_matrix %>% 
  select(c(rownames(tbspp_matrix))) %>% 
  # take out 0 sites
  mutate(total = rowSums(across(where(is.numeric)))) %>% 
  filter(total > 0) %>% 
  select(-total) %>% 
  as.matrix()

# making sure the rownames in the trait matrix == columns in the community matrix
rownames(tbspp_matrix) == colnames(cwm_commat)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
CWM_fixed <- functcomp(x = tbspp_matrix, a = cwm_commat, 
                  CWM.type = c("all")) %>% 
  rownames_to_column("sample_ID") %>% 
  rename_with(.fn = ~ paste0(., "_fixed", sep = ""), .cols = 2:28) %>% 
  left_join(., community_metadata, by = "sample_ID")

i. Thickness

# plot 
ggplot(CWM_fixed, aes(x = urchin_density, y = thickness_mm_fixed)) +
  geom_point(aes(col = site)) +
  geom_smooth(aes(col = site), method = "lm") +
  theme_bw()
## `geom_smooth()` using formula 'y ~ x'

ggplot(CWM_fixed, aes(x = thickness_mm_fixed)) +
  geom_histogram(binwidth = 0.05)

model_thickness_01a <- glmmTMB(thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + site + (1|year),
        data = CWM_fixed, family = gaussian(link = "identity")) 

model_thickness_01b <- glmmTMB(thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|year),
        data = CWM_fixed, family = gaussian(link = "identity")) 
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small eigenvalues detected. See
## vignette('troubleshooting')
model_thickness_01c <- glmmTMB(thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
        data = CWM_fixed, family = gaussian(link = "identity")) 

AIC(model_thickness_01a, model_thickness_01b, model_thickness_01c)
##                     df       AIC
## model_thickness_01a 11 -803.6142
## model_thickness_01b  7 -733.3296
## model_thickness_01c  8 -791.2190
# plot residuals
model_thickness_01a_resid <- simulateResiduals(model_thickness_01a)
plot(model_thickness_01a_resid)

model_thickness_01b_resid <- simulateResiduals(model_thickness_01b)
plot(model_thickness_01b_resid)

model_thickness_02_resid <- simulateResiduals(model_thickness_02)
plot(model_thickness_02_resid)

summary(model_thickness_01)
##  Family: gaussian  ( identity )
## Formula:          thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol +  
##     kelp_dry + site + (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   -800.2   -758.4    411.1   -822.2      321 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  year     (Intercept) 0.000482 0.02195 
##  Residual             0.004674 0.06837 
## Number of obs: 332, groups:  year, 14
## 
## Dispersion estimate for gaussian family (sigma^2): 0.00467 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      3.136e-01  1.571e-02  19.966  < 2e-16 ***
## urchin_density  -1.847e-03  5.515e-04  -3.349 0.000811 ***
## mean_sand        5.136e-04  3.068e-04   1.674 0.094111 .  
## mean_daily_umol -1.751e-06  9.921e-07  -1.765 0.077553 .  
## kelp_dry         2.230e-05  1.050e-05   2.123 0.033746 *  
## sitecarp         8.434e-02  1.485e-02   5.679 1.36e-08 ***
## siteivee         7.621e-03  1.457e-02   0.523 0.601022    
## sitemohk         6.433e-02  1.654e-02   3.889 0.000100 ***
## sitenapl         5.450e-02  1.459e-02   3.736 0.000187 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_thickness_02)
## Linear mixed model fit by REML ['lmerMod']
## Formula: thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol +  
##     kelp_dry + (1 | site) + (1 | year)
##    Data: CWM_fixed
## 
## REML criterion at convergence: -724.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1667 -0.7238 -0.0136  0.6721  3.3539 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  year     (Intercept) 0.0005743 0.02396 
##  site     (Intercept) 0.0012382 0.03519 
##  Residual             0.0047831 0.06916 
## Number of obs: 332, groups:  year, 14; site, 5
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)      3.546e-01  1.939e-02  18.292
## urchin_density  -1.727e-03  5.435e-04  -3.177
## mean_sand        4.464e-04  3.065e-04   1.457
## mean_daily_umol -1.770e-06  9.560e-07  -1.851
## kelp_dry         2.299e-05  1.061e-05   2.166
## 
## Correlation of Fixed Effects:
##             (Intr) urchn_ mn_snd mn_dl_
## urchn_dnsty -0.217                     
## mean_sand   -0.171  0.075              
## mean_dly_ml -0.235 -0.041  0.103       
## kelp_dry    -0.250  0.122  0.135 -0.146
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
# check collinearity
check_collinearity(model_thickness_01)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##             Term  VIF Increased SE Tolerance
##   urchin_density 1.20         1.10      0.83
##        mean_sand 1.56         1.25      0.64
##  mean_daily_umol 1.06         1.03      0.94
##         kelp_dry 1.64         1.28      0.61
##             site 2.36         1.54      0.42
plot(check_collinearity(model_thickness_01))

ii. LDMC/TDMC

ggplot(CWM_fixed, aes(x = urchin_density, y = tdmc_fixed)) +
  geom_point(aes(col = site)) +
  geom_smooth(aes(col = site), method = "lm")
## `geom_smooth()` using formula 'y ~ x'

ggplot(CWM_fixed, aes(x = tdmc_fixed)) +
  geom_histogram(binwidth = 0.01)

model_tdmc_01 <- glmmTMB(tdmc_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
        data = CWM_fixed, family = gaussian(link = "identity")) 

testDispersion(model_tdmc_01)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
## 
## data:  simulationOutput
## dispersion = 1.1352, p-value = 0.512
## alternative hypothesis: two.sided
testResiduals(model_tdmc_01)

## $uniformity
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.208, p-value = 1.335e-12
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
## 
## data:  simulationOutput
## dispersion = 1.1352, p-value = 0.512
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 2, observations = 324, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0007484357 0.0221193534
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                             0.00617284
## $uniformity
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.208, p-value = 1.335e-12
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
## 
## data:  simulationOutput
## dispersion = 1.1352, p-value = 0.512
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 2, observations = 324, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0007484357 0.0221193534
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                             0.00617284
model_tdmc_02 <- lmer(tdmc_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
        data = CWM_fixed)
## Warning: Some predictor variables are on very different scales: consider rescaling
# plot residuals
model_tdmc_01_resid <- simulateResiduals(model_tdmc_01)
plot(model_tdmc_01_resid)
## qu = 0.75, log(sigma) = -2.428833 : outer Newton did not converge fully.

model_tdmc_02_resid <- simulateResiduals(model_tdmc_02)
plot(model_tdmc_02_resid)
## qu = 0.75, log(sigma) = -1.775149 : outer Newton did not converge fully.

summary(model_tdmc_01)
##  Family: gaussian  ( identity )
## Formula:          tdmc_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry +  
##     (1 | site) + (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   -674.0   -643.8    345.0   -690.0      316 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.003701 0.06084 
##  year     (Intercept) 0.000248 0.01575 
##  Residual             0.006423 0.08015 
## Number of obs: 324, groups:  site, 5; year, 14
## 
## Dispersion estimate for gaussian family (sigma^2): 0.00642 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      2.498e-01  3.013e-02   8.292  < 2e-16 ***
## urchin_density  -1.221e-03  6.529e-04  -1.870   0.0615 .  
## mean_sand       -6.485e-04  3.582e-04  -1.811   0.0702 .  
## mean_daily_umol -5.158e-06  1.196e-06  -4.313 1.61e-05 ***
## kelp_dry        -1.119e-05  1.231e-05  -0.909   0.3631    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_tdmc_02)
## Linear mixed model fit by REML ['lmerMod']
## Formula: tdmc_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry +      (1 | site) + (1 | year)
##    Data: CWM_fixed
## 
## REML criterion at convergence: -611.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.85301 -0.55790 -0.05674  0.41836  2.68137 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  year     (Intercept) 0.000293 0.01712 
##  site     (Intercept) 0.004645 0.06816 
##  Residual             0.006486 0.08053 
## Number of obs: 324, groups:  year, 14; site, 5
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)      2.493e-01  3.261e-02   7.645
## urchin_density  -1.203e-03  6.296e-04  -1.910
## mean_sand       -6.431e-04  3.581e-04  -1.796
## mean_daily_umol -5.095e-06  1.036e-06  -4.919
## kelp_dry        -1.138e-05  1.240e-05  -0.918
## 
## Correlation of Fixed Effects:
##             (Intr) urchn_ mn_snd mn_dl_
## urchn_dnsty -0.133                     
## mean_sand   -0.117  0.051              
## mean_dly_ml -0.148 -0.122  0.131       
## kelp_dry    -0.184  0.136  0.125 -0.113
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

iii. SA:P ratio

ggplot(CWM_fixed, aes(x = urchin_density, y = ratio_fixed)) +
  geom_point(aes(col = site)) +
  geom_smooth(aes(col = site), method = "lm")
## `geom_smooth()` using formula 'y ~ x'

ggplot(CWM_fixed, aes(x = ratio_fixed)) +
  geom_histogram(binwidth = 1)

# full
model_ratio_01 <- glmmTMB(ratio_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
        data = CWM_fixed, family = gaussian(link = "identity")) 

# plot residuals
model_ratio_01_resid <- simulateResiduals(model_ratio_01)
plot(model_ratio_01_resid)

model_ratio_02_resid <- simulateResiduals(model_ratio_02)
plot(model_ratio_02_resid)

summary(model_ratio_01)
##  Family: gaussian  ( identity )
## Formula:          ratio_fixed ~ urchin_density + mean_sand + mean_daily_umol +  
##     kelp_dry + (1 | site) + (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   2016.8   2047.1  -1000.4   2000.8      316 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 16.084   4.011   
##  year     (Intercept)  3.826   1.956   
##  Residual             25.062   5.006   
## Number of obs: 324, groups:  site, 5; year, 14
## 
## Dispersion estimate for gaussian family (sigma^2): 25.1 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     7.209e+00  2.004e+00   3.597 0.000321 ***
## urchin_density  1.377e-01  4.130e-02   3.333 0.000858 ***
## mean_sand       9.882e-02  2.245e-02   4.401 1.08e-05 ***
## mean_daily_umol 3.901e-04  7.413e-05   5.262 1.42e-07 ***
## kelp_dry        8.832e-04  7.924e-04   1.115 0.265049    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_ratio_02)
## Linear mixed model fit by REML ['lmerMod']
## Formula: ratio_fixed ~ urchin_density + mean_sand + mean_daily_umol +  
##     kelp_dry + (1 | site) + (1 | year)
##    Data: CWM_fixed
## 
## REML criterion at convergence: 2061.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2541 -0.5628 -0.1867  0.5733  3.2969 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  year     (Intercept)  1.915   1.384   
##  site     (Intercept) 28.577   5.346   
##  Residual             25.241   5.024   
## Number of obs: 329, groups:  year, 14; site, 5
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)     6.9180720  2.5116439   2.754
## urchin_density  0.1225312  0.0397446   3.083
## mean_sand       0.1017156  0.0224186   4.537
## mean_daily_umol 0.0004240  0.0000678   6.253
## kelp_dry        0.0009980  0.0007681   1.299
## 
## Correlation of Fixed Effects:
##             (Intr) urchn_ mn_snd mn_dl_
## urchn_dnsty -0.112                     
## mean_sand   -0.096  0.053              
## mean_dly_ml -0.129 -0.086  0.124       
## kelp_dry    -0.146  0.122  0.131 -0.121
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
AIC(model_ratio_01, model_ratio_02)
## Warning in AIC.default(model_ratio_01, model_ratio_02): models are not all fitted to the same number of
## observations
##                df      AIC
## model_ratio_01  8 2016.829
## model_ratio_02  8 2077.731

iv. FvFm

ggplot(CWM_fixed, aes(x = urchin_density, y = fvfm_fixed)) +
  geom_point(aes(col = site)) +
  geom_smooth(aes(col = site), method = "lm")
## `geom_smooth()` using formula 'y ~ x'

ggplot(CWM_fixed, aes(x = fvfm_fixed)) +
  geom_histogram(binwidth = 0.01)

# full
model_fvfm_01a <- glmmTMB(fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
        data = CWM_fixed, family = Gamma(link = "logit")) 

# site as fixed instead of random effect
model_fvfm_01b <- glmmTMB(fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + site + (1|year),
        data = CWM_fixed, family = Gamma(link = "logit")) 

# no site
model_fvfm_01c <- glmmTMB(fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|year),
        data = CWM_fixed, family = Gamma(link = "logit")) 

# plot residuals
model_fvfm_01a_resid <- simulateResiduals(model_fvfm_01a)
plot(model_fvfm_01a_resid)

model_fvfm_01b_resid <- simulateResiduals(model_fvfm_01b)
plot(model_fvfm_01b_resid)

model_fvfm_01c_resid <- simulateResiduals(model_fvfm_01c)
plot(model_fvfm_01c_resid)

summary(model_fvfm_01a)
##  Family: Gamma  ( logit )
## Formula:          fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry +  
##     (1 | site) + (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   -698.9   -668.6    357.4   -714.9      316 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03093  0.1759  
##  year   (Intercept) 0.02150  0.1466  
## Number of obs: 324, groups:  site, 5; year, 14
## 
## Dispersion estimate for Gamma family (sigma^2): 0.0214 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -8.328e-02  9.826e-02  -0.848 0.396709    
## urchin_density   8.903e-03  2.766e-03   3.219 0.001287 ** 
## mean_sand        3.131e-03  1.368e-03   2.288 0.022137 *  
## mean_daily_umol  1.607e-05  4.741e-06   3.390 0.000699 ***
## kelp_dry         1.922e-05  4.719e-05   0.407 0.683823    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_fvfm_01b)
##  Family: Gamma  ( logit )
## Formula:          fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry +  
##     site + (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   -711.2   -669.7    366.6   -733.2      313 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  year   (Intercept) 0.01944  0.1394  
## Number of obs: 324, groups:  year, 14
## 
## Dispersion estimate for Gamma family (sigma^2): 0.0212 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -5.289e-02  7.433e-02  -0.712 0.476774    
## urchin_density   9.231e-03  2.757e-03   3.348 0.000813 ***
## mean_sand        3.127e-03  1.368e-03   2.285 0.022321 *  
## mean_daily_umol  1.666e-05  4.744e-06   3.512 0.000445 ***
## kelp_dry         2.211e-05  4.683e-05   0.472 0.636884    
## sitecarp        -2.969e-01  6.435e-02  -4.615 3.93e-06 ***
## siteivee        -4.220e-02  6.266e-02  -0.673 0.500657    
## sitemohk        -1.078e-01  7.375e-02  -1.462 0.143674    
## sitenapl         2.468e-01  6.493e-02   3.801 0.000144 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_fvfm_01c)
##  Family: Gamma  ( logit )
## Formula:          fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry +      (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   -618.8   -592.3    316.4   -632.8      317 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  year   (Intercept) 0.02717  0.1648  
## Number of obs: 324, groups:  year, 14
## 
## Dispersion estimate for Gamma family (sigma^2): 0.0288 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)  
## (Intercept)      5.947e-02  6.342e-02   0.938   0.3483  
## urchin_density   4.859e-03  3.094e-03   1.570   0.1163  
## mean_sand        2.123e-03  1.419e-03   1.496   0.1347  
## mean_daily_umol  8.661e-06  4.547e-06   1.905   0.0568 .
## kelp_dry        -4.695e-05  5.563e-05  -0.844   0.3987  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_fvfm_01a, model_fvfm_01b, model_fvfm_01c)
##                df       AIC
## model_fvfm_01a  8 -698.8904
## model_fvfm_01b 11 -711.2456
## model_fvfm_01c  7 -618.8032

v. max height

ggplot(CWM_fixed, aes(x = urchin_density, y = max_height_cm_fixed)) +
  geom_point(aes(col = site)) +
  geom_smooth(aes(col = site), method = "lm") +
  theme(legend.position = "none") +
  facet_wrap(~site)
## `geom_smooth()` using formula 'y ~ x'

ggplot(CWM_fixed, aes(x = max_height_cm_fixed)) +
  geom_histogram(binwidth = 1)

# full
model_max_height_cm_01a <- glmmTMB(max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
        data = CWM_fixed, family = gaussian(link = "identity")) 

# site as fixed instead of random effect
model_max_height_cm_01b <- glmmTMB(max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + site + (1|year),
        data = CWM_fixed, family = gaussian(link = "identity")) 

# no site
model_max_height_cm_01c <- glmmTMB(max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|year),
        data = CWM_fixed, family = gaussian(link = "identity")) 

# plot residuals
model_max_height_cm_01a_resid <- simulateResiduals(model_max_height_cm_01a)
plot(model_max_height_cm_01a_resid)

model_max_height_cm_01b_resid <- simulateResiduals(model_max_height_cm_01b)
plot(model_max_height_cm_01b_resid)

model_max_height_cm_01c_resid <- simulateResiduals(model_max_height_cm_01c)
plot(model_max_height_cm_01c_resid)

summary(model_max_height_cm_01a)
##  Family: gaussian  ( identity )
## Formula:          max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol +  
##     kelp_dry + (1 | site) + (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   2348.2   2378.5  -1166.1   2332.2      316 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 30.43    5.517   
##  year     (Intercept) 12.02    3.467   
##  Residual             69.83    8.357   
## Number of obs: 324, groups:  site, 5; year, 14
## 
## Dispersion estimate for gaussian family (sigma^2): 69.8 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     2.045e+01  2.899e+00   7.055 1.73e-12 ***
## urchin_density  2.239e-01  6.910e-02   3.240   0.0012 ** 
## mean_sand       1.558e-01  3.744e-02   4.162 3.16e-05 ***
## mean_daily_umol 5.403e-04  1.243e-04   4.347 1.38e-05 ***
## kelp_dry        7.651e-04  1.324e-03   0.578   0.5635    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_max_height_cm_01b)
##  Family: gaussian  ( identity )
## Formula:          max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol +  
##     kelp_dry + site + (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   2334.4   2376.0  -1156.2   2312.4      313 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  year     (Intercept) 10.87    3.297   
##  Residual             68.97    8.305   
## Number of obs: 324, groups:  year, 14
## 
## Dispersion estimate for gaussian family (sigma^2):   69 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     16.9548387  1.9815632   8.556  < 2e-16 ***
## urchin_density   0.2277658  0.0687214   3.314 0.000919 ***
## mean_sand        0.1585724  0.0373990   4.240 2.24e-05 ***
## mean_daily_umol  0.0005753  0.0001225   4.695 2.66e-06 ***
## kelp_dry         0.0008155  0.0013166   0.619 0.535630    
## sitecarp        -2.6336951  1.8341814  -1.436 0.151032    
## siteivee         6.5766697  1.8002934   3.653 0.000259 ***
## sitemohk        -0.5999541  2.0136300  -0.298 0.765744    
## sitenapl        12.5192031  1.7950750   6.974 3.08e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_max_height_cm_01c)
##  Family: gaussian  ( identity )
## Formula:          max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol +  
##     kelp_dry + (1 | year)
## Data: CWM_fixed
## 
##      AIC      BIC   logLik deviance df.resid 
##   2432.7   2459.1  -1209.3   2418.7      317 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  year     (Intercept) 19.52    4.419   
##  Residual             94.88    9.741   
## Number of obs: 324, groups:  year, 14
## 
## Dispersion estimate for gaussian family (sigma^2): 94.9 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      2.571e+01  1.696e+00  15.160  < 2e-16 ***
## urchin_density   1.504e-01  7.831e-02   1.920  0.05485 .  
## mean_sand        1.011e-01  3.872e-02   2.610  0.00904 ** 
## mean_daily_umol -8.313e-06  1.205e-04  -0.069  0.94501    
## kelp_dry        -2.262e-04  1.502e-03  -0.151  0.88025    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_max_height_cm_01a, model_max_height_cm_01b, model_max_height_cm_01c)
##                         df      AIC
## model_max_height_cm_01a  8 2348.228
## model_max_height_cm_01b 11 2334.373
## model_max_height_cm_01c  7 2432.673

b. specific CWM using FD::dbFD()

i. Bullito

Only one species, skip.

ii. Arroyo Quemado

Only one species, skip.

iii. Isla Vista

# get a site specific trait by species matrix
ivee_tbspp <- specif_tbspp("Isla Vista")

# get a site specific transect by species matrix
ivee_commat <- specif_commat("ivee", ivee_tbspp)

# making sure the rownames in the trait matrix == columns in the community matrix
rownames(ivee_tbspp) == colnames(ivee_commat)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
# calculate CWM for site
CWM_ivee <- functcomp(x = ivee_tbspp, a = ivee_commat, 
                  CWM.type = c("all")) %>% 
  rownames_to_column("sample_ID") %>% 
  rename_with(.fn = ~ paste0(., "_specif", sep = ""), .cols = 2:ncol(.)) %>% 
  left_join(., CWM_fixed, by = "sample_ID") %>% 
  mutate(site = "ivee")

iv. Mohawk

# get a site specific trait by species matrix
mohk_tbspp <- specif_tbspp("Mohawk")

# get a site specific transect by species matrix
mohk_commat <- specif_commat("mohk", mohk_tbspp)

# making sure the rownames in the trait matrix == columns in the community matrix
rownames(mohk_tbspp) == colnames(mohk_commat)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# calculate CWM for site
CWM_mohk <- functcomp(x = mohk_tbspp, a = mohk_commat, 
                  CWM.type = c("all")) %>% 
  rownames_to_column("sample_ID") %>% 
  rename_with(.fn = ~ paste0(., "_specif", sep = ""), .cols = 2:ncol(.)) %>% 
  left_join(., CWM_fixed, by = "sample_ID") %>% 
  mutate(site = "mohk")

v. Carpinteria

Only one species, skip.

c. putting all the data frames together

df <- CWM_mohk %>% 
  select(sample_ID, thickness_mm_specif, tdmc_specif, site) %>% 
  rbind(., (CWM_ivee %>% select(sample_ID, thickness_mm_specif, tdmc_specif, site))) %>% 
  rbind(., (CWM_bull %>% select(sample_ID, thickness_mm_specif, tdmc_specif, site))) %>% 
  left_join(., CWM_fixed, by = "sample_ID")

d. variance decomp using cati::traitflex.anova()

tdmc_tf <- traitflex.anova(~ site.x, specif.avg = tdmc_specif, const.avg = tdmc_fixed, data = df)
tdmc_tf
## 
##  Decomposing trait sum of squares into composition turnover
##  effect, intraspecific trait variability, and their covariation:
##             Turnover Intraspec. Covariation  Total
## site.x      0.091056     0.1393     0.11144 0.3418
## Residuals   0.772176     1.4119    -0.52404 1.6600
## Total       0.863232     1.5512    -0.41260 2.0018
## 
##  Relative contributions:
##             Turnover Intraspec. Covariation  Total
## site.x       0.04549    0.06959     0.05567 0.1707
## Residuals    0.38574    0.70530    -0.26178 0.8293
## Total        0.43122    0.77489    -0.20611 1.0000
## 
##  Significance of testable effects:
##               Turnover Intraspec.     Total
## site.x      0.00016746 0.00064961 4.141e-07
plot(tdmc_tf)

thickness_tf <- traitflex.anova(~ site.x, specif.avg = thickness_mm_specif, const.avg = thickness_mm_fixed, data = df)
thickness_tf
## 
##  Decomposing trait sum of squares into composition turnover
##  effect, intraspecific trait variability, and their covariation:
##             Turnover Intraspec. Covariation   Total
## site.x       0.12055    0.26870     0.16142 0.55067
## Residuals    0.56746    0.55187    -0.63954 0.47979
## Total        0.68801    0.82056    -0.47812 1.03046
## 
##  Relative contributions:
##             Turnover Intraspec. Covariation  Total
## site.x        0.1170     0.2608      0.1566 0.5344
## Residuals     0.5507     0.5356     -0.6206 0.4656
## Total         0.6677     0.7963     -0.4640 1.0000
## 
##  Significance of testable effects:
##               Turnover Intraspec.      Total
## site.x      2.9824e-07 3.6505e-14 8.7061e-27
plot(thickness_tf)